An inverse-free ADI algorithm for computing Lagrangian invariant subspaces

The numerical computation of Lagrangian invariant subspaces of large scale Hamiltonian matrices is discussed in the context of the solution of Lyapunov and Riccati equations. A new version of the low-rank alternating direction implicit method is introduced, which in order to avoid numerical difficulties with solutions that are of very large norm, uses an inverse-free representation of the subspace and avoids inverses of ill-conditioned matrices. It is shown that this prevents large growth of the elements of the solution which may destroy a low-rank approximation of the solution. A partial error analysis is presented and the behavior of the method is demonstrated via several numerical examples.